Οικολογική και Στατιστική Μοντελοποίηση

1 Εγκατάσταση βιβλιοθήκης

Εγκατάσταση του πρώτου και κύριου πακέτου (και ενός βοηθητικού)
(Αυτό γίνεται μια φορά σε κάθε Η/Υ - εάν έχει εγκατασταθεί δεν χρειάζεται να δοθεί ξανά η εντολή αυτή)

## CRAN repository:
install.packages("pacman", repos = "http://cran.us.r-project.org")

2 Φόρτωση βιβλιοθήκης

Φόρτωση του εν λόγω πακέτου

library(pacman)

3 Εγκατάσταση βιβλιοθηκών

Με την εντολή p_load('insert_package_name_here') γίνεται η εγκατάσταση των πακέτων που μας ενδιαφέρουν

Δοκιμάστε το

4 Φόρτωση βιβλιοθηκών

Εντολή για φόρτωση των πακέτων

pacman::p_load("betapart", "maptree", "ape", "ggdendro", "dendextend", "plyr", 
    "zoo", "ggplot2", "vegan", "magrittr", "gdm", "devtools", "tidyr", "dplyr", 
    "picante")

devtools::install_github("daijiang/rtrees")

source("https://bioconductor.org/biocLite.R")

biocLite("ggtree")

## if (!requireNamespace('BiocManager', quietly = TRUE))
## install.packages('BiocManager')

## BiocManager::install('ggtree')

5 Working Directory

Με την ακόλουθη εντολή, βλέπουμε σε ποιόν φάκελο δουλεύουμε:

getwd()

6 Αλλαγή του working directory

Ενώ με την εντολή αυτή, αλλάζουμε τον φάκελο στον οποίο δουλεύουμε και αποθηκεύουμε δεδομένα
[π.χ., setwd("E:/") εάν θέλουμε να δουλεύουμε στον σκληρό δίσκο Ε]

setwd("E:/Academic/Lessons/EMB/Labs/Labs")

7 Εισαγωγή αρχείου

Τώρα ας εισάγουμε το αρχείο το οποίο περιέχει τα δεδομένα τα οποία μας ενδιαφέρουν να αναλύσουμε

plants <- read.csv("G:/Academic/Lessons/EMB/Labs/Labs/Ryukyus.csv", h = T, sep = ";", 
    row.names = 1)
dim(plants)

Όπως θα διαπιστώσετε, σήμερα θα ασχοληθούμε πάλι με τα Ryukyus Islands, τα οποία είναι περισσότερο γνωστά ως Nansei Islands (στα Ιαπωνικά: Nansei-shoto, δηλαδή “Νοτιοδυτικά Νησιά”), ανήκουν στην Ιαπωνία. Αποτελούν ένα νησιωτικό τόξο, το οποίο βρίσκεται σε γεωγραφικό πλάτος 24ο - 31.0οΝ και εκτείνεται σε μήκος 1300 km σε έναν άξονα ΝΔ-ΒΑ, ξεκινώντας από το Kyushu και καταλήγει στην Taiwan. Το νοτιότερο νησί είναι το Yonaguni, ενώ το μεγαλύτερο η Okinawa. Τα μεγαλύτερα σε μέγεθος νησιά είναι κυρίως ορεινά, ενώ τα μικρότερα είναι κοραλλιογενή. Καθώς εντοπίζεται στην ζώνη μετάβασης από το τροπικό στο θερμό εύκρατο κλίμα (δηλαδή, βρίσκεται στην υποτροπική ζώνη), οι κλιματικές συνθήκες είναι ήπιες στο αρχιπέλαγος αυτό, καθ’όλη την διάρκεια του έτους, με μέση θερμοκρασία 15οC το χειμώνα και 28οC το καλοκαίρι. Βέβαια, υπάρχουν διαβαθμίσεις: στο Β τμήμα του αρχιπελάγους, το κλίμα χαρακτηρίζεται ως υγρό υποτροπικό, ενώ στο Ν τμήμα του ως τροπικό (ταξινόμηση κατά Koppen). Εξαιτίας της γεωγραφικής τους θέσης, στα νησιά αυτά δεν υπάρχει ξηρή περίοδος και τα περισσότερα εξ αυτών δέχονται περισσότερα από 2000 m βροχής, λόγω των μουσώνων και των καλοκαιρινών τυφώνων. Ως εκ τούτου, στα Ryukyus Islands δεν εμφανίζεται η υποτροπική ξηρή ζώνη, η οποία διαχωρίζει τα εύκρατα δάση από τις υγρές τροπικές περιοχές σε άλλα μέρη του πλανήτη και η φυσική τους βλάστηση αποτελείται από πλατύφυλλα αειθαλή δαση (Good, 1845; Maekawa, 1974; Kira, 1991).

Γεωλογικά, το αρχιπέλαγος διαιρείται σε 3 τμήματα:
1. Northern Ryukyus (τα νησιά βορείως του Tokara Strait)
2. Central Ryukyus (τα νησιά μεταξύ του Tokara Strait και του Kerama Gap)
3. Southern Ryukyus (τα νησιά νοτίως του Kerama Gap)

Εικόνα 1. Γεωτεκτονική διαίρεση του αρχιπελάγους Ryukyus.

Εικόνα 1. Γεωτεκτονική διαίρεση του αρχιπελάγους Ryukyus.

Τα Ryukyus είναι ηπειρωτικά νησιά και αποτελούν τις κορυφές της κορδιλλιέρας του αρχιπελάγους Ryukyus και έχουν δημιουργηθεί ως αποτέλεσμα της καταβύθισης της τεκτονικής πλάκας των Φιλιππίνων (Ota, 1998).

8 Κατασκευή μήτρας ανομοιότητας

Για τη διερεύνηση των παραγόντων που επηρεάζουν τα πρότυπα της β-ποικιλότητας, θα χρειαστεί να μορφοποιήσουμε κατάλληλα τα δεδομένα τα οποία έχουμε στην διάθεση μας.

Στο σημείο αυτό θέλουμε να δούμε πόσο μοιάζουν (ή όχι) οι κοινότητες των νησιών μεταξύ τους. Για τον σκοπό αυτό, θα δημιουργήσουμε μια μήτρα ανομοιότητας. Υπάρχουν πολλοί τρόποι (π.χ., μέσω της βιβλιοθήκης vegan ή της βιβλιοθήκης betapart) και πολλοί δείκτες (δείκτης Sorensen ή δείκτης Jaccard). Θα βασιστούμε στην β-ποικιλότητα (εναλλαγή ειδών.). Σήμερα, θα κατασκευάσουμε μια μήτρα ανομοιότητας βάσει της β-ποικιλότητας και του δείκτη Jaccard, τόσο για την ταξινομική, όσο και για την φυλογενετική όψη της ποικιλότητας των ειδών.

Προς το παρόν, ας κατασκευάσουμε μια μήτρα ανομοιότητας:

## ============================================================================##
## Create a taxonomic dissimilarity matrix
## ============================================================================##

beta.plants <- beta.pair(plants[, 1:1815], index.family = "jaccard")

## beta.jtu --> Turnover

## beta.jne --> Nestedness

## beta.jac --> Total beta diversity

str(beta.plants)
## List of 3
##  $ beta.jtu: 'dist' num [1:325] 0.5302 0.2283 0.0595 0.0361 0.3452 ...
##   ..- attr(*, "Labels")= chr [1:26] "Takeshima" "Ioujima" "Kuroshima" "Tanegashima" ...
##   ..- attr(*, "Size")= int 26
##   ..- attr(*, "call")= language as.dist.default(m = beta.jtu)
##   ..- attr(*, "Diag")= logi FALSE
##   ..- attr(*, "Upper")= logi FALSE
##  $ beta.jne: 'dist' num [1:325] 0.0107 0.4541 0.7533 0.8082 0.2831 ...
##   ..- attr(*, "Labels")= chr [1:26] "Takeshima" "Ioujima" "Kuroshima" "Tanegashima" ...
##   ..- attr(*, "Size")= int 26
##   ..- attr(*, "call")= language as.dist.default(m = beta.jne)
##   ..- attr(*, "Diag")= logi FALSE
##   ..- attr(*, "Upper")= logi FALSE
##  $ beta.jac: 'dist' num [1:325] 0.541 0.682 0.813 0.844 0.628 ...
##   ..- attr(*, "Labels")= chr [1:26] "Takeshima" "Ioujima" "Kuroshima" "Tanegashima" ...
##   ..- attr(*, "Size")= int 26
##   ..- attr(*, "call")= language as.dist.default(m = beta.jac)
##   ..- attr(*, "Diag")= logi FALSE
##   ..- attr(*, "Upper")= logi FALSE

9 Κατασκευή φυλογενετικού δένδρου

Προκειμένου να μπορέσουμε να υπολογίσουμε την φυλογενετική β-ποικιλότητα για το αρχιπέλαγος της Ιαπωνίας, θα χρειαστεί να κατασκευάσουμε πρώτα το φυλογενετικό δένδρο των ειδών που απαρτίζουν τις νησιωτικές φυτοκοινότητες. Για τον σκοπό αυτόν, θα χρησιμοποιήσουμε εντολές από τις βιβλιοθήκες tidyr, dplyr & rtrees. Η τελευταία βιβλιοθήκη περιλαμβάνει τα φυλογενετικά δένδρα για αρκετούς φυτικούς οργανισμούς, καθώς και για τους ιχθύες, τα πτηνά και τα θηλαστικά. Αρχικά, θα μετατρέψουμε την μήτρα παρουσίας/απουσίας σε έναν πίνακα που θα περιέχει το όνομα του είδους και του γένους για κάθε φυτικό taxon το οποίο απαντάται στα Ryukyus και στη συνέχεια, θα χρησιμοποιήσουμε αυτόν τον πίνακα για να κατασκευάσουμε το φυλογενετικό μας δένδρο. Ακόμα και εάν το φυλογενετικό μας δένδρο περιέχει πολυτομίες, είναι ασφαλές να το χρησιμοποιήσουμε για μακροοικολογικές αναλύσεις (Li et al., 2019).

##============================================================================##
## Create the phylogenetic tree
##============================================================================##
library(rtrees)
jpn_tree <- plants %>% 
  pivot_longer(Codonacanthus_pauciflorus:Alpinia_japonica, 
               names_to = 'species') %>% 
  dplyr::select(species) %>% 
  distinct(species) %>% 
  separate(species, 'genus', 
           sep = '_', 
           remove = F) %>% 
  get_tree(sp_list = .,
           tree = tree_plant_otl, # either 
           taxon = "plant", # or
           scenario = "S1",
           show_grafted = TRUE)
##============================================================================##

Καθώς κάποια από τα ονόματα των ειδών (tip.labels) στο φυλογενετικό δένδρο που μόλις κατασκευάσαμε περιέχουν το σύμβολο *, θα χρησιμοποιήσουμε την ακόλουθη εντολή για να αφαιρέσουμε το εν λόγω σύμβολο:

## ============================================================================##
## Delete the asterisk
## ============================================================================##
jpn_tree$tip.label <- gsub("\\*", "", jpn_tree$tip.label)
## ============================================================================##

Στη συνέχεια, μπορούμε να αναπαραστήσουμε γραφικά το φυλογενετικό δένδρο, δίνοντας διαφορετικό χρώμα σε κάθε γένος:

## ============================================================================##
## Plot the tree
## ============================================================================##
library(ggtree)
groupInfo <- split(jpn_tree$tip.label, gsub("_\\w+", "", jpn_tree$tip.label))

groups_jpn <- groupOTU(jpn_tree, groupInfo)
ggtree(groups_jpn, aes(color = group), layout = "circular") + geom_tiplab2(size = 0.5, 
    align = T)

9.1 Υπολογισμός φυλογενετικής β-ποικιλότητας

Στο σημείο αυτό, θα χρησιμοποιήσουμε μια εντολή από την βιβλιοθήκη picante, προκειμένου να κρατήσουμε μόνο εκείνα τα είδη τα οποία είναι κοινά μεταξύ του φυλογενετικού δένδρου και της μήτρας παρουσίας/απουσίας:

## ============================================================================##
## Combine the trees with the P/A data
## ============================================================================##
combined.japan <- picante::match.phylo.comm(jpn_tree, plants[, 1:1815])  #check for mismatches/missing species


## the resulting object is a list with $phy and $comm elements replace our
## original data with the sorted/matched data
phy.japan <- combined.japan$phy

comm.japan <- combined.japan$comm
## ============================================================================##

Τώρα μπορούμε να υπολογίσουμε την φυλογενετική β-ποικιλότητα, μέσω της βιβλιοθήκης betapart:

## ============================================================================##
## Compute phylogenetic beta diversity
## ============================================================================##
pbd <- phylo.beta.pair(comm.japan, phy.japan, index.family = "jaccard")
pbd_turnover <- pbd$phylo.beta.jtu %>% as.matrix()
pbd_nestedness <- pbd$phylo.beta.jne %>% as.matrix()
pbd_overall <- pbd$phylo.beta.jac %>% as.matrix()
## ============================================================================##

9.2 Μεταφόρτωση αβιοτικών δεδομένων και συντελεστής συσχέτισης

Αρχικά θα φορτώσουμε τα αβιοτικά μας δεδομένα και εν συνεχεία θα ελέγξουμε ποιά εξ αυτών εμφανίζουν υψηλό συντελεστή συσχέτισης, προκειμένου να τα αφαιρέσουμε από την ανάλυση:

## ============================================================================##
## Load the data
## ============================================================================##
jap_isl <- readxl::read_excel("G:/Academic/Lessons/EMB/Labs/Labs/Ryukyu islands/Ryukyus.xlsx")
## ============================================================================##


## ============================================================================##
## Check the correlation structure
## ============================================================================##
usdm::vifcor(jap_isl %>% as.data.frame() %>% dplyr::select(-c(1:3)), th = 0.7)
## ============================================================================##

Στη συνέχεια, μπορούμε να φορτώσουμε τα αβιοτικά δεδομένα τα οποία δεν εμφανίζουν πρόβλημα συσχέτισης:

## ============================================================================##
## Load the final IVs
## ============================================================================##
envTab <- readxl::read_excel("G:/Academic/Lessons/EMB/Labs/Labs/Ryukyu islands/Ryukyus 2.xlsx") %>% 
    as.data.frame()
## ============================================================================##

Ακολούθως, μπορούμε να συνδυάσουμε την μήτρα ανομοιότητας που δημιουργήσαμε προηγουμένως, τόσο για την ταξινομική, όσο και για την φυλογενετική β-ποικιλότητα, με τα αβιοτικά μας δεδομένα, ώστε να δημιουργήσουμε - με την απαραίτητη μορφή για την ανάλυση - τα εξαρτημένα μας δεδομένα:

## ============================================================================##
## Rename the beta diversity dissimilarity matrix
## ============================================================================##
beta_div <- beta.plants$beta.jac %>% as.matrix()

beta_div_phylo <- pbd_overall

## ============================================================================##


## ============================================================================##
## Join and create the DV
## ============================================================================##
sppTab <- envTab %>% dplyr::select(Longtitude:Site) %>% full_join(., beta_div %>% 
    as.data.frame() %>% mutate(Site = envTab$Site))

sppTab_phylo <- envTab %>% dplyr::select(Longtitude:Site) %>% full_join(., pbd_overall %>% 
    as.data.frame() %>% mutate(Site = envTab$Site))

## ============================================================================##


## ============================================================================##
## Remove the first two columns
## ============================================================================##
sppTab <- sppTab[, -c(1, 2)]
sppTab_phylo <- sppTab_phylo[, -c(1, 2)]
## ============================================================================##

Πλέον είμαστε σε θέση να αναλύσουμε τα δεδομένα μας:

## ============================================================================##
## Run the analysis
## ============================================================================##
gdmTab <- formatsitepair(bioData = sppTab, bioFormat = 3, siteColumn = "Site", 
    XColumn = "Longtitude", YColumn = "Latitude", predData = envTab, sampleSites = 1)

gdmModel.1 <- gdm(gdmTab, geo = T, splines = NULL, knots = NULL)


gdmTab_phylo <- formatsitepair(bioData = sppTab_phylo, bioFormat = 3, siteColumn = "Site", 
    XColumn = "Longtitude", YColumn = "Latitude", predData = envTab, sampleSites = 1)

gdmModel.1_phylo <- gdm(gdmTab_phylo, geo = T, splines = NULL, knots = NULL)

## ============================================================================##

Ας δούμε το ποσοστό της ποικιλότητας που εξηγείται από τις ανεξάρτητες μεταβλητές μας:

## ============================================================================##
## Deviance explained
## ============================================================================##
str(gdmModel.1)
## List of 15
##  $ dataname    : symbol gdmTab
##  $ geo         : logi TRUE
##  $ sample      : int 300
##  $ gdmdeviance : num 21.6
##  $ nulldeviance: num 26.1
##  $ explained   : num 17.1
##  $ intercept   : num 0.7
##  $ predictors  : chr [1:7] "Geographic" "Area" "bio_3" "bio_5" ...
##  $ coefficients: num [1:21] 0 0 0.0903 0.1916 0.2245 ...
##  $ knots       : num [1:21] 0.132 2.191 9.501 0.965 34.398 ...
##  $ splines     : num [1:7] 3 3 3 3 3 3 3
##  $ creationdate: chr "Sat May 23 18:59:09 2020"
##  $ observed    : num [1:300] 0.807 0.564 0.569 0.786 0.75 ...
##  $ predicted   : num [1:300] 0.706 0.597 0.574 0.712 0.714 ...
##  $ ecological  : num [1:300] 1.223 0.909 0.853 1.244 1.253 ...
##  - attr(*, "class")= chr [1:2] "gdm" "list"
str(gdmModel.1_phylo)
## List of 15
##  $ dataname    : symbol gdmTab_phylo
##  $ geo         : logi TRUE
##  $ sample      : int 300
##  $ gdmdeviance : num 19.7
##  $ nulldeviance: num 23.1
##  $ explained   : num 14.7
##  $ intercept   : num 0.507
##  $ predictors  : chr [1:7] "Geographic" "Area" "bio_3" "bio_5" ...
##  $ coefficients: num [1:21] 0 0 0.0534 0.1126 0.1559 ...
##  $ knots       : num [1:21] 0.132 2.191 9.501 0.965 34.398 ...
##  $ splines     : num [1:7] 3 3 3 3 3 3 3
##  $ creationdate: chr "Sat May 23 18:59:11 2020"
##  $ observed    : num [1:300] 0.676 0.438 0.452 0.639 0.604 ...
##  $ predicted   : num [1:300] 0.573 0.476 0.457 0.574 0.575 ...
##  $ ecological  : num [1:300] 0.852 0.647 0.61 0.853 0.856 ...
##  - attr(*, "class")= chr [1:2] "gdm" "list"

Μπορούμε επίσης να δούμε με ποιόν τρόπο επηρεάζει η κάθε μεταβλητή την β-ποικιλότητα:

## ============================================================================##
## Plot the results
## ============================================================================##
plot(gdmModel.1, plot = c(2, 3), plot.color = "blue", plot.linewidth = 2, include.rug = TRUE)

plot(gdmModel.1_phylo, plot = c(2, 3), plot.color = "blue", plot.linewidth = 2, 
    include.rug = TRUE)

Μπορούμε να διαπιστώσουμε ποια ομάδα μεταβλητών (π.χ., κλιματικές, γεωγραφικές ή άλλες μεταβλητές) επηρεάζει σε μεγαλύτερο βαθμό τα πρότυπα της β-ποικιλότητας:

## ============================================================================##
## Source the function
## ============================================================================##
source("G:/Academic/Lessons/EMB/Labs/Labs/GDM - Scotland/Scripts/Scripts from Matt/gdm.partition.deviance.R")
## ============================================================================##


## ============================================================================##
## First make the varSet argument
## ============================================================================##
varSet <- vector("list", 2)  ## makes a list of length 2
## ============================================================================##


## ============================================================================##
## Names of climate variables
## ============================================================================##
varSet[[1]] <- c("bio_5", "bio_8", "bio_14", "bio_18")
## ============================================================================##


## ============================================================================##
## names of topography variables
## ============================================================================##
varSet[[2]] <- "Area"
## ============================================================================##


## ============================================================================##
## Now add names to the list to indicate that variable types
## ============================================================================##
names(varSet) <- c("climate", "area")
## ============================================================================##


## ============================================================================##
## Check the result
## ============================================================================##
(scPart <- gdm.partition.deviance(sitePairTable = gdmTab, varSets = varSet, 
    partSpace = T))
VARIABLE_SET DEVIANCE
3 area 8.688
2 climate 12.415
1 geo 3.226
4 climate and area 16.992
5 geo and climate 12.475
6 geo and area 10.642
7 ALL VARIABLES (geo and climate and area) 17.122
8 area alone 4.647
9 climate alone 6.480
10 geo alone 0.130
11 climate union area, exclude geo 2.769
12 climate union geo, exclude area 1.825
13 geo union area, exclude climate -0.070
14 area union climate union geo 1.342
15 Unexplained 82.878 
(scPart <- gdm.partition.deviance(sitePairTable = gdmTab_phylo, varSets = varSet, 
    partSpace = T))
VARIABLE_SET DEVIANCE
3 area 7.580
2 climate 11.094
1 geo 2.061
4 climate and area 14.606
5 geo and climate 11.124
6 geo and area 8.754
7 ALL VARIABLES (geo and climate and area) 14.697
8 area alone 3.572
9 climate alone 5.942
10 geo alone 0.091
11 climate union area, exclude geo 3.121
12 climate union geo, exclude area 1.084
13 geo union area, exclude climate -0.061
14 area union climate union geo 0.948
15 Unexplained 85.303

Τέλος, μπορούμε να διαπιστώσουμε ποια μεταβλητή είναι η σημαντικότερη όσον αφορά την διαμόρφωση των προτύπων της β-ποικιλότητας:

##============================================================================##
## Variable importance
##============================================================================##
gdmModel.1.sig <- gdm.varImp(gdmTab,
                             geo = T,
                             fullModelOnly = FALSE,
                             nPerm = 999, ## Use 999 always
                             parallel = T)

gdmModel.1.sig_phylo <- gdm.varImp(gdmTab_phylo,
                             geo = T,
                             fullModelOnly = FALSE,
                             nPerm = 999, ## Use 999 always 
                             parallel = T)

##============================================================================##

Πρώτα θα δούμε τη σημαντικότητα του μοντέλου μας:

## ============================================================================##
## Disable scientific notation
## ============================================================================##
options(scipen = 999)
## ============================================================================##


## ============================================================================##
## Model significance
## ============================================================================##
gdmModel.1.sig[[1]]
##                              fullModel  fullModel-1  fullModel-2
## Model deviance              21.6102271  21.61022714  21.66356117
## Percent deviance explained  17.1218956  17.12189557  16.91735247
## Model p-value                0.1681682   0.09859155   0.09264854
## Fitted permutations        999.0000000 994.00000000 993.00000000
##                             fullModel-3  fullModel-4  fullModel-5
## Model deviance              21.71214545  21.89205698  23.09391255
## Percent deviance explained  16.73102527  16.04104055  11.43176407
## Model p-value                0.09090909   0.09176225   0.06626506
## Fitted permutations        990.00000000 959.00000000 830.00000000
##                             fullModel-6
## Model deviance              23.42164842
## Percent deviance explained  10.17485326
## Model p-value                0.07492795
## Fitted permutations        694.00000000
gdmModel.1.sig_phylo[[1]]
##                              fullModel fullModel-1 fullModel-2 fullModel-3
## Model deviance              19.6685399  19.6685399  19.7008017  19.7690448
## Percent deviance explained  14.6965921  14.6965921  14.5566712  14.2606977
## Model p-value                0.2112112   0.1533066   0.1395582   0.1123482
## Fitted permutations        999.0000000 998.0000000 996.0000000 988.0000000
##                            fullModel-4 fullModel-5 fullModel-6
## Model deviance              19.8963951   20.725409  20.9462932
## Percent deviance explained  13.7083732   10.112900   9.1549144
## Model p-value                0.1129707    0.113879   0.1022727
## Fitted permutations        956.0000000  843.000000 704.0000000

Στη συνέχεια θα δούμε τη σημαντικότητα των μεταβλητών μας:

## ============================================================================##
## Variable importance
## ============================================================================##
gdmModel.1.sig[[2]]
##             fullModel fullModel-1 fullModel-2 fullModel-3 fullModel-4
## Geographic  0.7565835   0.7565835   0.9379181    1.746937    4.781302
## Area       27.1412341  27.1412341  26.2925957   26.459044   28.734274
## bio_3       0.0000000          NA          NA          NA          NA
## bio_5       1.1946288   1.1946288          NA          NA          NA
## bio_8       2.7261250   2.7261250   4.0474748    4.123984          NA
## bio_14      1.4387203   1.4387203   1.1013969          NA          NA
## bio_18     25.6612032  25.6612040  29.6398932   29.475639   33.657629
##            fullModel-5
## Geographic    10.99490
## Area                NA
## bio_3               NA
## bio_5               NA
## bio_8               NA
## bio_14              NA
## bio_18        71.77729
gdmModel.1.sig_phylo[[2]]
##             fullModel fullModel-1 fullModel-2 fullModel-3 fullModel-4
## Geographic  0.6188968   0.6188968    0.751926    1.783489    4.477375
## Area       24.3073448  24.3073448   23.580249   24.068295   26.228298
## bio_3       0.0000000          NA          NA          NA          NA
## bio_5       0.9520638   0.9520638          NA          NA          NA
## bio_8       2.7157448   2.7157448    3.804445    3.873054          NA
## bio_14      2.3876666   2.3876666    2.033250          NA          NA
## bio_18     24.8346840  24.8346840   31.244601   31.170221   36.140487
##            fullModel-5
## Geographic    9.472909
## Area                NA
## bio_3               NA
## bio_5               NA
## bio_8               NA
## bio_14              NA
## bio_18       79.618075

Ας την αναπαραστήσουμε γραφικά:

## ============================================================================##
## Plot them
## ============================================================================##
barplot(sort(gdmModel.1.sig[[2]][, 1], decreasing = T))

barplot(sort(gdmModel.1.sig_phylo[[2]][, 1], decreasing = T))

Τέλος, μπορούμε να δούμε το p-value της κάθε μεταβλητής:

##============================================================================##
## Each variable's p-value
##============================================================================##
gdmModel.1.sig[[3]]
gdmModel.1.sig_phylo[[3]]

Βιβλιογραφία

Good, R. (1845) The geography of the flowering plants. Longsmans,

Kira, T. (1991) Forest ecosystems of east and southeast asia in a global perspective. Ecological Research, 6, 185–200.

Li, D., Trotta, L., Marx, H.E., Allen, J.M., Sun, M., Soltis, D.E., Soltis, P.S., Guralnick, R.P., & Baiser, B. (2019) For common community phylogenetic analyses, go ahead and use synthesis phylogenies. Ecology, 100, e02788.

Maekawa, F. (1974) Origin and characteristics of japan’s flora. The Flora and Vegetation of Japan., 33–86.

Ota, H. (1998) Geographic patterns of endemism and speciation in amphibians and reptiles of the ryukyu archipelago, japan, with special reference to their paleogeographical implications. Researches on Population Ecology, 40, 189–204.

Κώστας Κουγιουμουτζής

Spring semester 2019-2020

 

Delivered to you by Dr. Kostas Kougioumoutzis

kkougiou@aua.gr